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O ■ ABSTRACT 



The temperature and ionization of SN 1987A is modeled between 200 and 2000 
days in its nebular phase, using a time-dependent model. We include all important 
Q ■ elements, as well as the primary composition zones in the supernova. The energy 

r — . input is provided by radioactive decay of 56 Co, 57 Co, and 44 Ti. The thermalization 

of the resulting gamma-rays and positrons is calculated by solving the Spencer- Fano 
equation. Both the ionization and the individual level populations are calculated 
time-dependently Adiabatic cooling is included in the energy equation. Charge 
£vi I transfer is important for determining the ionization, and is included with available and 

estimated rates. Full, multilevel atoms are used for the observationally important ions. 
As input models to the calculations we use explosion models for SN 1987A calculated 



C . by Woosley et al. and Nomoto et al. 

The most important result in this paper refers to the evolution of the temperature 



and ionization of the various abundance zones. The metal-rich core undergoes a 
thermal instability, often referred to as the IR-catastrophe, at 600 — 1000 days. The 



o . . . 

hydrogen-rich zones evolve adiabatically after 500 — 800 days, while in the helium 
region both adiabatic cooling and line cooling are of equal importance after ~ 1000 
days. Freeze-out of the recombination is important in the hydrogen and helium zones. 
Concomitant with the IR-catastrophe, the bulk of the emission shifts from optical and 
near-IR lines to the mid- and far-IR. After the IR-catastrophe, the cooling is mainly 
due to far-IR lines and adiabatic expansion. 

Dust cooling is likely to be important in the zones where dust forms. We find that 
the dust condensation temperatures occur later than ~ 500 days in the oxygen-rich 
zones, and the most favorable zone for dust condensation is the iron core. The 
uncertainties introduced by the, in some cases, unknown charge transfer rates are 
discussed. Especially for ions of low abundances differences can be substantial. 



Subject headings: abundances - line formation - nucleosynthesis - supernovae: general 
- supernovae: individual (SN 1987A) - stars: evolution - stars: interiors 



1. INTRODUCTION 

The observations of SN 1987A represent a unique set of data. Especially at late epochs 
high quality spectra exist at epochs when other supernovae are too faint to be observed even 
photometrically. This is particularly true for the IR region (e.g., Chevalier 1992, McCray 1993). 
The analysis of these observations are, however, not trivial and represent a considerably more 
complex situation compared to that of a planetary nebula, or an H II region. Complications 
include large optical depths of even forbidden lines, high velocities and consequent severe line 
blending, non-thermal excitation and ionization, non- uniform abundances, molecule and dust 
formation, as well as time-dependent effects. For reviews of the physics see McCray (1993) and 
Fransson (1994). In this paper and in Kozma & Fransson (1997) (Paper II) we concentrate on the 
spectral modeling in the nebular phase when the core region is transparent in at least the optical 
and IR regions. In a subsequent paper we will discuss the broad-band photometry, including the 
bolometric light curve. 

Previous discussions of the late, nebular spectrum can be divided into attempts of self- 
consistent modeling of the spectrum, or analyses of the lines of particular elements or molecules. 
Models in the first category are hampered in many cases by lack of atomic data, as well as the 
general complexity of the problem. Examples of this type of models include some early predictions 
by Fransson & Chevalier (1987) and Colgan & Hollenbach (1988), and modeling of observed 
spectra by Swartz, Harkness, & Wheeler (1989). Applications to core collapse supernovae in 
general, including Type Ib's and Type Ic's can be found in Fransson & Chevalier (1989) and 
Swartz (1994). The 'diagnostic' approach has been discussed in a series of papers by Li & McCray 
(1992, 1993, 1995), Li, McCray & Sunyaev (1993), Xu et al. (1992), Chugai (1994), and by Kozma 
& Fransson (1992, hereafter KF92). While the problem in this case is simplified to an easily 
controllable problem, it is often limited by the lack of knowledge of temperature, ionization and 
radiation field. As will be apparent in this paper and in Paper II, the limitation to one specific 
chemical composition, and in most cases a single element or ion, is especially severe. 

Predictions are one of the best methods for judging the success of a complex model. It is 
therefore of some interest to discuss the paper by Fransson & Chevalier (1987), where immediately 
after the explosion predictions for the spectral evolution were made. There are several limitations 
to that paper. The most serious were the neglect of the hydrogen zone, and the simplified 
treatment of the iron emission. The model was therefore only applicable to the helium, oxygen 
and silicon - sulphur regions. Even with these limitations there were several interesting results. 
An important prediction in this paper was the dramatic drop in temperature at ~ 600 — 700 
days in the metal-rich regions, and consequent transition from optical to far-IR radiation. This is 
commonly referred to as the IR-catastrophe (e.g., Axelrod 1980, Fransson 1994). The optical to 
far-IR transition is one of the characteristic features in the evolution of SN 1987A, and is clearly 
reflected if the light curves of the lines are plotted as a fraction of the total bolometric luminosity 
(e.g., Menzies 1991). Often this has been attributed to dust thermalization, which is certainly 
important, but independent of this, just the normal cooling in the metal-rich core will have this 



effect. The drop in temperature is also most likely a requirement for molecule and dust formation, 
and was in fact anticipated by Fransson & Chevalier. The concomitant drop in the optical lines 
and the dust formation, as seen in the line profiles (Lucy et al. 1989, 1991) and mid-IR emission 
(e.g., Wooden et al. 1993), are therefore most likely connected. 

The predictions by Fransson & Chevalier have been directly compared to the observations of 
the light curves of the strongest lines, by Menzies (1991), Spyromilio et al. (1991), and Meikle et al. 
(1993). While there are quantitative discrepancies, it was found that the general shape of e.g. the 
strong [Ca II] A A 7291, 7324 and [0 1] A A 6300,6364 light curves agreed well qualitatively with 
the predictions. Note that most of the [Ca II] emission in this model comes from the oxygen-rich 
zone, and not from the calcium-rich Si - S zone, showing that a small amount of an effective 
cooler, like Ca II, can dominate the emission from a zone with a large gamma-ray deposition. In 
Paper II we will see that this is also the case in these models, although it is the hydrogen-rich 
zones, rather than the oxygen-rich wich dominate the emission, in agreement with Li & McCray 
(1993). 

Another interesting result in the paper by Fransson & Chevalier was the non-uniform 
temperature and ionization in the different composition zones, with the metal-rich regions having 
considerably lower temperature. This illustrates the necessity to take all important abundance 
components into account, as well as the importance of distinguishing between macroscopic mixing 
of blobs of different compositions, compared to microscopic mixing of the gas. 

An important aspect of the late evolution of the spectrum and light curve is the departure 
from (quasi-) steady state. As density decreases oc t -3 , recombination and cooling time scales 
increase. Eventually, as these become longer than the radioactive decay time scale, r, the spectrum 
and light curve have to be calculated time-dependently. In addition, when the radiative cooling 
time becomes longer than the expansion time scale, adiabatic cooling must be included. These 
effects were discussed in Fransson &; Kozma (1993), with particular emphasis on the consequences 
for the bolometric light curve. Implicit in this calculation was, however, a complete calculation 
of the spectrum as a function of time. Some results of this were reported in Fransson, Houck & 
Kozma (1996), and it was shown especially that the [0 1] A A 6300, 6364 lines provided compelling 
evidence for the IR-catastrophe. In addition, the importance of non-thermal excitation for the 
spectrum at late time were stressed. These points are discussed further in Paper II. 

In these two papers we continue our attempt to model the evolution of the line emission 
from SN 1987A as realistically as possible. Consequently, we calculate the temperature and 
ionization self-consistently, using a time-dependent formalism. Non-thermal processes are included 
with the method in KF92. In addition, we include all important abundance zones in the ejecta 
with appropriate filling factors. The important ions are modeled as multilevel atoms, using the 
best available atomic data. Even with these ingredients there are a number of deficiencies in 
our modeling. Probably the most serious omission is that we do not include the resonance line 
blanketing, especially important in the UV. Further, the molecular chemistry and cooling are 



not included. Although we have tried to include the best atomic data available, we are still 
handicapped by a lack of reliable data, especially with regard to charge transfer reactions and 
photoionization of the iron group elements. In Paper II we argue that these omissions, although 
important for specific lines, should not influence our general conclusions severely. We hope to 
improve upon these shortcomings in the future. 

In this paper we discuss our model and the temperature and ionization evolution in detail. In 
Paper II we concentrate on the line emission, as well as on the line profiles. In section g we discuss 
our basic model, assumptions, and physical processes included, as well as the explosion models we 
use. Section |3| contains a detailed discussion of our results with regard to the temperature and 
ionization, while in section || we discuss the sensitivity of our results to especially the UV-radiation 
field and the poorly known charge transfer rates. In section || we discuss the dust formation and 
the effects of dust cooling. In section || our main results are for convenience summarized. In the 
appendix we have collected the references for the atomic data. 



2. THE MODEL 

In these papers we focus on the emission from SN 1987A at times later than ~ 200 days, 
when the supernova has entered the nebular phase. As the matter is expanding and becoming 
more transparent, it is possible to extract information on the nucleosynthesis, which has taken 
place in the progenitor and in the explosion, as well as on the hydrodynamics of the explosion. 

Our model calculates the temperature, ionization and level populations of the different regions 
in the supernova envelope. The first time step is calculated in steady state, which is a good 
approximation at 200 days. The following time steps are calculated time-dependently. Knowing 
the time evolution of the physical properties of the supernova matter we calculate the bolometric 
luminosity, as well as luminosities for ~ 5900 lines. As input, the calculations require an explosion 
model, as well as the radial structure of the ejecta with abundances and masses as a function of 
velocity. Our model is spherical symmetric, with zones of different compositions and thicknesses. 
Mixing in velocity is taken into account in an approximate way by mixing shells from different 
abundance zones in velocity. 

We have the option to switch between either doing the calculations in steady state or 
time-dependently. The steady state calculations are done primarily in order to check that energy 
is conserved at every time step. We find that energy is conserved to within ~ 1 % at all times. 
By comparing the steady state and time-dependent calculations it is possible to determine the 
level of freeze-out. The temperature, ionization, and level populations of the multi-level atoms are 
calculated time-dependently. 

Below, we first discuss the physical processes included in our model, and then the explosion 
models used in our calculations. 



2.1. Physical Processes 

2.1.1. Temperature Calculation 

The temperature, T, of a region is determined by balancing the cooling and heating 
mechanisms. In steady state the total cooling (Ann e ) at a certain point is just balanced by the 
total heating (H), 

K{T)nn e = H. (1) 

From this condition it is possible to determine the temperature. In the time-dependent case the 
temperature change can be expressed as 

-(H - A(T)nn e ) — 



alt 3k(l + x e )n t 1 + x e dt 

where n is the number density of ions and atoms, n e the number density of free electrons, and 
x e = n e /n. The second term in this equation describes adiabatic cooling. If this term dominates, 

T oc t- 2 . 

Most of the cooling is by collisional excitation of bound-bound transitions. The total line 
cooling is a sum over all the transitions i,j, 

L hncs = ^2(n i n e Ci j (T)-n j n e C j i(T))hiy ij , (3) 

where Cij(T) and Cji(T) are the collisional excitation, and de-excitation rates, respectively 
(Osterbrock 1989), ni is the number densities of ions in level i, and hi/ij is the photon energy. Line 
cooling is the dominant cooling process, but we also include cooling due to recombinations and 
free-free cooling. We approximate recombination cooling with 

Lrcc ~ Yl nk+in e a k kT , (4) 

k 

where the summation is over all ions k, and at refers to the recombination from k + 1 to k. 
Free-free cooling is given by 

L ff = Y, L426 xW~ 27 T 1/2 Z k l n e n k <g s >. (5) 

k 

The summation is over the ions k, Z k is the charge of ion k, and < gs > is the mean gaunt factor 
for free-free emission. 

Cooling, both due to lines, recombinations, and free-free emission, is in the low density, 
optically thin limit, proportional to the electron fraction. We therefore define a cooling function, 
A(T), so that the total cooling equals A(T)nn e = -Lii ncs + L rcc + Lg, 

HT) = SEtxiC^T) - XjC ii {T))hv i j+ 

k ij 



^x k+1 a k kT + 1.426 x KT^T 1 / 2 ^ < g s >] , (6) 

where x% is the fraction of ions in level i, and the summation over k represents all ions. A depends 
sensitively on the electron temperature, T. 

The heating in the ejecta is dominated by non-thermal heating, i? 7 , which is independent of 
the electron temperature. This process is discussed in detail in KF92. We express the non-thermal 
heating as 

i? 7 = 4-7T Jya^rjan , (7) 

where \tt J 7 <7 7 is the deposited non-thermal energy per particle. This quantity is discussed in 
more detail below, tjh is the fraction of the deposited energy going into heating the free, thermal 
electrons. ??h depends on the electron fraction, and increases with increasing x e (KF92). 

Photoelectric heating from the ground states is included for all ions k, and from excited levels 
i, where there is atomic data available. The total photoelectric heating can be expressed, 

^p = EE n kti r ^ J ^-^ du , (8) 

k i Juk > i hv 

where hv k ^ is the ionization energy, a v is the photoionization cross section from level % in ion A;, 
J u is the mean radiation intensity at u, and Xk,i is the ionization potential. 

The total heating is therefore 

H = H y + H P . (9) 



2.1.2. Ionization Balance 

In our calculations we include H I-II, He I-III, C I-III, N III, O I-III, Ne I-II, Na III, Mg 
I-III, Si I-III, S I-III, Ar I-II, Ca I-III, Fe I-V, Co II, and Ni I-II. 

For an estimated value of the electron fraction, x e , the ionization balance for an element with 
only two ionization stages, k and k + 1 is solved from 

-^T + — = -{^k + Pk)n k -n k 2^niE, kl + 

a k n e n k+1 + n k+1 ^ n^fc+l,«- ( 10 ) 

i 

For elements with three or more ionization stages equation ( |l0| ) is generalized to include processes 
both below and above the relevant ionization stage. Here, T k is the non-thermal ionization rate, 

Fk = 47T Ji^L , (ii) 

Xeff,k 



where Xeff,k is the effective ionization potential, defined as Xeff,k = XkXk/Vk (see KF92), \k being 
the ionization potential, x& the number fraction, and % the fraction of the deposited non-thermal 
energy going into ionization of ion k. P^ is the photoionization rate, 

p '= 4 '£^" • < 12 > 

while ctfc is the recombination coefficient from k + 1 to k. In our calculations Pk and atk include 
photoionizations from, and recombinations to excited levels in ion k. 

In many cases charge transfer reactions with ions I may contribute, either as an ionization, 
n kni£,kJ ') or recombination, rik+iniQJTn, process (Fransson & Chevalier 1989, Swartz 1994). 
Excited levels of ion k may also be involved. 

Finally, the number conservation equation is 

/2 nk = "element- (13) 

k 

Equations (|10|) and (|13| ) are solved using the Newton- Raphson method. For stability and efficiency, 
equations ( |10|) and (|13| ) are solved implicitly with backward differencing, 



■n 



n+l 



At 



/W 1 ) ( 14 ) 



where n£ + refers to the new value and n£ to the old value. The function f(n% + ) on the right 
hand side includes all ionization and recombination processes in equation ([[(])• Even though 
the implicit treatment requires iteration, it makes the overall problem more stable, and allows 
considerably longer time steps. The time step varies between the different zones, depending 
on the cooling and recombination time scales. The hydrogen and helium-rich regions therefore 
require smaller number of time steps, compared to the metal-rich regions, where these time scales 
are short. For the same reason, the late stages allow longer time steps, compared to the first 
epochs. As an example, the iron zones require ~ 15000 time steps from 200 - 2000 days, while the 
corresponding numbers are ~ 4000 for the oxygen zones, and ~ 1000 for the hydrogen zones. 



2.1.3. Level Populations 

With the abundance of each ionization stage known, the level populations are calculated 
for a number of ions. For the following multilevel atoms we calculate the level populations 
time-dependently : HI (30 levels with n < 20), Hel (16 levels), OI (13 levels), Call (6 levels), 
Fel (121 levels), Fell (191 levels), Fe III (110 levels), and Fe IV (43 levels). Lines from other 
elements are solved, either as two-level atoms, or, in the case of forbidden lines, as 4-6 level 
systems (see appendix). 



The rate equations for the level populations are given by 

dn k ,i . 3n k>i 

— Jj- H ^— = -rik,i{-rk,i + J- k,i + Oii,k-irie) 

~ n k,i / j \fHjPij ~r ^ij n e + *ij) 
j<i 

+ Yl n k,j(,Cjin e + Tji + Pji) 

j<i 

j>i 

+ ^^ nkj(Ajipji + Gjin e + /jj) 
+n fe+ ia fc ,jn e + n fe _i(P i)fe _i + r iife _i) 

i i 

+n k ^J2^k-i,i (15) 

i 

where the index i refers to the actual level, and j to other levels within the ion, k to the ionization 
stage of the element, and I to any other element. P kt i and T k ,i are the photoionization rate and 
non-thermal ionization rate, from ion k, level i. o^fc— 1 is the recombination coefficient from ion k, 
and level i to ion k — 1. Downward and upward radiative, and collisional rates, Aij/3ij and Cijn e , 
are included. Here /3jj is the escape probability, which includes both the Sobolev escape probability 
and a continuum destruction probability. This is discussed in more detail in sections 2.1.5 and 
2~Q . 

Similar to the non-thermal ionization rate (eq. p| ) the non-thermal excitation rate T^ from 
the lower level i to the upper level j is given by 

r, - 47rJ T a T j (16) 

Xij x i/Vij 

where Xy is the excitation energy, rr, the number fraction in level i, and rjij is the fraction of the 
deposited energy going into excitation of the transition ij. For those elements where non-thermal 
excitations are included, only excitations from the ground state, i = 1, are calculated except for O 
I, Fel and Fell. For these, non-thermal excitations from the individual levels of the ground state 
multiplet are included. 

Pij and Pji are the photoexcitation/deexcitation rates between levels i and j. The Einstein 
relations and the Sobolev assumption (e.g., Tielens & Hollenbach 1985) result in 

c 2 
= 2hv 3 



and 



Pij — — Pji- (18) 

9i 



In accordance to our neglect of line blanketing, we only include the continuum contribution to 
Jji (see section 2.1. 7| ) . Further, a^j is the recombination coefficient from ion k + 1 to an individual 



level i of ion k. For ions once or more ionized, photoionization and non-thermal ionization from 
ion k — 1 are included. We assume that these ionizations end up in the first level, i = 1, Pj=i &-i 
and Fj = i fe_i, with the exception of O II. For oxygen we include non-thermal ionization from O I 
to the excited states 2 D and 2 P , as well as to the ground state of O II. 

Charge transfer reactions with other elements may either act as an ionization or recombination 
process. We assume that charge transfer reactions occur from their ground states, i = 1, with the 
exception for He I. Swartz (1994) gives rates for charge transfer involving excited levels in He I. 
The term n^n^J represents the charge transfer rate for the reaction A k + B l — > A k±l + B lzfl , and 
n k±i n iCk±i i t ne charge transfer rate for the reaction A + B > — > A + B . Charge transfer 



reactions are further discussed in section 2.1 



Together with the particle conservation equation 

re fc = H n 'M' ( 19 ) 

i 

the rate equations are again solved implicitly, using Newton-Raphson iteration. 

2.1.4- Non-Thermal Deposition 

The energy source for the ejecta at late times is well established to be the decay of radioactive 
isotopes created in the explosion. The dominant isotope synthesized is 56 Ni. 56 Ni rapidly decays 
to 56 Co on a time scale of Tdecay = 8.8 days. 56 Co then decays to 56 Fe on a longer time scale of 
111.26 days. In the decay of 56 Co 96.5 % of the energy is emitted in the form of gamma-rays while 
the remaining 3.5 % of the energy is emitted as positrons. Other important radioactive isotopes 
created in the explosion are 57 Ni and 44 Ti. 57 Ni decays rapidly to 57 Co, and 57 Co to 57 Fe with 
Tdecay = 391.2 days. Decays of 57 Co only result in gamma-rays, no positrons. 44 Ti decays first to 
44 Sc on a time scale of 83.7 years followed by a rapid decay to 44 Ca. In the decay of 44 Sc positrons 
are also emitted. 

The gamma-ray luminosity due to the decay of 56 Co is 



; M( 56 Ni) ^ /m . 26d 
0.1 M Q 



L 7 = 1.26 X 1Q 42 ^V ^V e -t/111.26« erg g -l (20) 



of 57 Co, 

L 1 = 1.37 x 10 41 M( 57 Ni) e - 4 / 391 - 2d erg s" 1 , (21) 



10 



and of 44 Ti, 



L 7 = 4.1 x 10 36 M{ ^T) e- f / 83 - 7y " erg s" 1 . (22) 



10" 4 M m 



These expressions include the gamma-rays due to positron annihilation. 

Positrons are assumed to be deposited locally in the iron-rich regions. This is supported even 
at 7.8 years by the analysis of the spectrum of SN 1987A by Chugai et al. (1997). The luminosity 
due to the energy input from positrons from 56 Co decay is 

^^"^S^ ergs " (23) 



and from Ti, 



L e+ = 1.3 x 10 36 M( r?? e-*/ 83 - 7Kr erg s" 1 . (24) 

10~ 4 M B v ; 



The life time of 44 Ti is uncertain, with values ranging from 66.9 - 96.1 years (see Timmes et 
al. 1996). Here we use the latest value of 83.7 ± 15 years (Meissner et al. 1995). In our calculations 
we include 0.07 M of 56 Ni, from the bolometric light curve, 3.3 x 10 -3 M of 57 Ni, which 
corresponds to a 57 Ni/ 56 Ni ratio equal to 2 times the solar 57 Fe/ 56 Fe ratio, and 10 -4 M Q of 44 Ti 
(Woosley, Pinto & Weaver 1988, Kumagai et al. 1991, Woosley & Hoffman 1991). These masses 
correspond to our standard case in Fransson & Kozma (1993). With these masses, 56 Co dominates 
the energy input up to ~ 1200 days. Thereafter, 57 Co takes over up to ~ 1900 days, after which 
44 Ti dominates. 22 Na with a decay time of 3.75 years never contributes more than 10 % before 
2000 days (Woosley, Pinto & Hartman 1989). 

Thermalization of the gamma-rays and positrons is discussed in detail in KF92 and Li, 
McCray, & Sunyaev (1993). Using the Spencer- Fano formalism (e.g., Spencer & Fano 1954, 
Douthat 1975a,b, Xu 1989) we calculate the fractions of the deposited energy going into heating, 
ionization, and excitation, respectively. These fractions depend on the composition of the gas, the 
state of ionization and the electron fraction. The Spencer-Fano calculations are repeated for each 
composition zone, when the electron fraction has changed by more than 5 %, or by more than 0.01 
in absolute value. 

The gamma-ray intensity depends on the distribution of the radioactive elements. This is 
uncertain, and depends on the degree of mixing in the gas (see section |2.2j ). We distribute the 
radioactive source within the core, proportional to the iron mass. The ionization, excitation, and 
heating depend on 47rJ 7 <7 7 (erg s" 1 particle" 1 ), the total gamma-ray energy deposited per second 
per particle at a certain radius in the supernova ejecta. With particles we refer to ions and atoms 
(not electrons). J 7 (r) is the gamma-ray mean intensity at radius r, c 7 (r) (cm 2 particle -1 ) is the 
energy- averaged cross section for gamma-ray absorption at radius r for an ion. For the gamma-rays 
resulting from the decay of 56 Co, Fransson & Chevalier (1989) find k 7 = 0.06 (Zi/Ai) cm 2 g , 
where Zi and A{ is the nuclear charge and the atomic weight, respectively, for ion i. For 
hydrogen-rich regions Zi/Ai ~ 1, while in the other regions, dominated by helium or heavier 
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elements, Zi/Ai ~ 0.5. Therefore, we use k 7 ( 56 Co) = 0.06 cm 2 g _1 in the hydrogen-rich regions, 
and k 7 ( 56 Co) = 0.03 cm 2 g _1 in the remaining regions. The average cross section for a particular 
composition is then calculated from cr 7 ( 56 Co) = j4 mea n m p k 7 ( 56 Co), where j4 mea n is the mean 
atomic weight in the region. 

The mean energy of the gamma-rays emitted in the 57 Co-decay is lower than in the decay of 
56 Co. Woosley et al. (1989) find that k 7 ( 57 Co) as 2.4 k 7 ( 56 Co), and also an effective opacity for 
the gamma-rays from 44 Ti, K 7 ( 44 Ti) = 0.040, compared to k 7 ( 56 Co) = 0.033. In our calculations 
we use the the same opacity for 44 Ti and 56 Co. 

In the case of a central energy source the total deposited energy may be calculated from, 

4vra 7 J 7 (r) = [(1 - e- A ^ 56(r) )e- T ^ 56(r) L 7 , 56 + 

(1 - e- Ar ^ 57(r) )e- T ^ 57(r ' ) L 7i 5 7 + 

(1 _ e -Ar 7 , 44 M )e -r 7 , 44 (r) L ^ 44]/n . on(r)) 

(25) 

where the index 56 refers to 56 Co, 57 to 57 Co, and 44 to Ti, the quantity r(r) is given below in 
equation (J27|), Ar(r) = a(r)n(r)Ar is the optical depth of the shell of thickness Ar, and nion(^) is 
the number of particles (i.e. ions) at r. 

If the radioactive isotopes are distributed throughout the core region, the calculation of J 7 
is more complicated. Knowing the luminosities for the radioactive isotopes (eqs. ]2(| - [22]), and 
their distribution (from the explosion model) we find the emissivity as a function of radius, j^(r) 
We then calculate the total deposited energy from 

47ra 7 J 7 (r) = 4vr- / / [^ f M r )hM^~ TlM) 

2 J -l Jo 

+a lM {r)j lM (l)e- T -<< i4{l) ] dl dfi. (26) 

For each r, inside the core (r < R CO re), we integrate over all angles, [i = cos 9, integration interval 
[-1,1]. For each [i we integrate along a path, I, from r (I = 0) to the border of the core region, i? C orc 
(I = 'max)- Outside of the core (r > R CO re) the integration interval for /i is [-1,-V^ 2 — Rcore/r]. 
The optical depth, r, in equation ( |2"6| ) is given by 

T (l) = f a(l')n(l')dl', (27) 

Jo 

for each isotope. 



In general we can write 
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^"leAL' (28) 



where i? CO re is the core radius. Limiting forms of D 7 are discussed in KF92 



2.1.5. Line Transfer 

In our calculations we take the radiative transfer of the lines into account by using the 
Sobolev approximation (Sobolev 1957, 1960, Castor 1970). The supernova ejecta is expanding 
homologously at the times we are interested in. Therefore, each mass element, with mass 
coordinate M, has a constant velocity, v(M) = constant, or r{M) = v(M)t, where t is the 
time since core collapse. The expansion velocity is much larger than the thermal velocity of the 
matter, and the Sobolev approximation is therefore valid for an individual, well separated line. 
With this approximation the line transfer becomes a purely local phenomena, either the photon 
is re-absorbed locally or escapes the medium. A photon emitted in a certain transition may, 
however, be absorbed by a line with a longer wavelength at a different point in the ejecta. This 
line blocking is important in the UV, and we hope to treat this in a subsequent paper (see also 
Fransson 1994, Li & McCray 1996, Houck & Fransson 1996). The effect of this has been checked 
by comparing calculations for SN1993J (Houck & Fransson 1996) with and without these processes 
included, where it was found that the strongest lines from dominant ions with high ionization 
potentials, like O I, were only slightly affected, while elements with low ionization potentials and 
low abundances, like Na I, are sensitive to the UV-field. The electron density and temperature 
were only marginally affected. We test the sensitivity of our results to the UV-field by having the 
option of switching it on or off. This is discussed in Paper II. Our crude treatment of the radiative 
transfer is one of the major weaknesses of our calculations. 

The Sobolev escape probability, /3 S , can be written 

1 - e~ Ts 
& = — . (29) 

T s 

Here r s is the Sobolev optical depth for a particular transition, ij (i < j), 

Ts = A 1 fii>§ L g lt{l _mn L) 

87T Qi gj m 

In addition of being re-absorbed in the line itself, a line photon may also be destroyed by 
continuum absorption, i.e., photoionization. Note that the Sobolev expression for the escape 
probability (eq. |30t| ) is independent of assumptions of the line profile (e.g., Rybicki 1984). 
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2.1.6. Continuum Destruction 

The continuum destruction probability, j3 c , has in the Sobolev approximation been calculated 
by Hummer & Rybicki (1985) for the Doppler case and complete redistribution, 

In this equation kc/kj^, is the ratio between the continuum opacity, kc, and the line core opacity, 
k^, and £ = l/r s . For a Doppler profile Hummer & Rybicki (1985) approximate F(£, ^) by 

F(£, g) « 2{- In 0r(|| + 2£[- In 0r(^ + 2e* )] 1 /2) } i/2, (32) 

where xo ~ 2 is the frequency displacement from the line center in thermal Doppler widths, 
(y — v)/ /\v-q. The assumptions made for this estimate are that kc/ki, <C 1 and £ <C 1. 

For strong resonance lines the damping wings may become important. In this case the 
continuum destruction probability above is modified. The destruction probability in the Voigt 
case with complete redistribution has for a static geometry been calculated by Hummer (1968), 
giving expressions and tabulations for the function F above. For kc/k^ ^ 0.1a and a <; 0.1, where 
a is the damping parameter, Hummer gives the approximation 



For fcc/^L ^ 0.1a equation (32) gives an adequate approximation. 



The validity of complete redistribution versus partial redistribution has been discussed 
extensively for solar conditions in the case of Lya, Ca II H & K and Mg II (e.g., Mihalas 1978, and 
references therein). It is there concluded that at densities lower than those prevailing in the outer 
chromosphere, partial redistribution is a better approximation to the line profile in the damping 
wings of the lines. Because the densities are in our case considerably lower than this, it is likely 
that partial redistribution is a more realistic assumption for the strong resonance lines from the 
ejecta. The continuum destruction has for partial redistribution been discussed by Chugai (1987), 
who finds the approximation 

This expression takes only scattering in the wings into account, and does not apply to the Doppler 
core, where complete redistribution, as given by equation (|32|), is adequate. 



The most important lines for which partial redistribution applies are Lya, He I A 584, Ca II 
H and K, O I A 1302, and Mg II A 2800, as well as the Fe II resonance lines. As will be seen in 
Paper II, the choice of destruction probability has important effects for especially the He I A 584 
line and the branching to the A 2.058 /im line. The effects of this are discussed in Paper II. 
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In our calculations we use the effective transition probability, Aji(3ji where the escape 
probability (3ji includes both the Sobolev escape probability, and the continuum destruction 
probability, 

&» = (!- &)& + &■ (35) 



2.1.7. Photoionization and Photoexcitation/deexcitation 

In our calculation of the ionization balance we include photoionization of all the ions with 
ionization potential smaller than 30 eV. Where data is available, we also include photoionization 
from excited levels in our multi-level atoms. Lines, recombination emission, two-photon emission 
and free-free emission may contribute to the photoionization. 

We assume that only the continuum emission (recombination, two-photon, and free-free 
emission) contribute to the photoexcitation/deexcitation. Line photons destroyed by continuum 
absorption (see previous section) are added to the photoionization of the different elements, 
according to their fractional opacity. References for the photoionization cross sections are given in 
the appendix. 



2.1.8. Charge Transfer and Penning Ionization 

Charge transfer reactions may be important in determining the ionization stages for the 
elements (Fransson 1994, Swartz 1994). This is especially the case for trace elements. The 
ionization balance of the dominating element in a zone is, however, not affected much by charge 
transfer reactions. Therefore, the total electron fraction is rather insensitive to these reactions. 
References to the sources of data are given in the appendix. 

In modeling the helium emission from SN 1987A the Penning ionization process, 

Hel (2 3 S) + HI -> Hel (l 1 ^) + HII + e~ (36) 

may be important for de-populating the excited level 2 3 S* of Hel (Chugai 1991, Li & McCray 
1995). We include this process in our model. 



2.2. The Explosion Model 

As input to our calculations we require a model for the composition and mass (or equivalently, 
the density) as a function of velocity. The assumption of homologous expansion relates the radius 
of a mass element to the velocity, r{M) = v(M)t. Explosion models for SN 1987A have been 
calculated by Nomoto & Hashimoto (1988), Woosley (1988), Arnett (1987), and Hashimoto, 
Nomoto, &: Shigeyama (1989). These models are all one-dimensional calculations. However, 
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several groups have found that the core is unstable to hydro dynamical instabilities, causing 
mixing of the various burning shells, as well as modifying the density distribution (e.g., Fryxell, 
Miiller, & Arnett 1991, Hachisu et al. 1992, Herant & Benz 1991, 1992). This mixing should be 
macroscopic, with little mixing on the microscopic scale (Fryxell et al. 1991). As emphasized by 
e.g. Fransson & Chevalier (1989) and Li, McCray & Syunyaev (1993), macroscopic mixing, rather 
than microscopic, can have a major influence on the observed spectrum. 

In our calculations we assume an overall spherical symmetry, with shells of different 
compositions and masses. Unfortunately, most models which treat the nucleosynthesis in sufficient 
detail are completely one-dimensional, with shells containing the heavy elements close to the 
center, and with continuously decreasing mean atomic weight as one proceeds outward in radius. 
Conversely, multi-dimensional models start with artificial initial perturbations, and do not treat 
the nucleosynthesis in detail. Because of limited computing power it is also not yet realistic to 
include multi-dimensional effects in the calculation of the spectrum, especially with regard to the 
radiative transfer. 

For these reasons we mimic the macroscopic mixing by using the results of the one-dimensional 
explosion models with regard to the nucleosynthesis, but not in density and velocity. Instead, 
we divide the ejecta into a core region, containing the metal-rich region, and where mixing is 
determining the distribution in velocity, and an envelope, containing most of the hydrogen, as 
well as a major part of the helium, and where mixing is less important. To simulate the structure 
of the multi-dimensional models in the core we divide this into a number of shells with varying 
compositions, placed in an alternating order. The density, pi, of each of these shells is determined 
by the total mass of the component in the core, Mi, and its filling factor, /j, so that 

_ 3 Mj 

A -4vr/ 4 y c 3 ore t3- W 

Each component is then divided into a number of shells, usually three, all having this density, 
placed uniformly within the core. In some cases (see below) we determine the location, and masses 
of these shells, so that a best fit to the line profiles is obtained. A schematic representation of the 
structure, including only a few of the zones, is given in Figure |]. 

Based on line widths of [O I], Mg I], and [Fe II] lines we take the velocity of the core to 
be 2000 km s" 1 . The structure of the envelope is based on the line profiles of Ha and the He 
I A 2.058 /urn line (Paper II). In our calculations we use a maximum expansion velocity for the 
envelope of 6000 km s . 

To represent the most important burning zones in the post-explosion ejecta we include six 
different abundance components. The dominant elements (by number) are Fe - He, Si - S, O - 
Si - S/O - Ne - Mg, O - C, He, and H. Abundances are either from the 10H model (Woosley 
& Weaver 1986, Woosley 1988) or the 11E1 model (Nomoto, & Hashimoto 1988, Shigeyama, 
Nomoto, & Hashimoto 1988, Shigeyama & Nomoto 1990, Hashimoto et al. 1989). The abundances 
used, together with the masses of the different components, for the 10H and 11E1 models are 
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presented in Tables [l] and ^. The sodium abundance is not given in the 10H model. We estimate 
this by scaling it from the magnesium abundance, so that the integrated X( 23 Na)/A"( 24 Mg) ratio 
is that given by the 11E1 model. This gives X(Na) ~ 8.2 x 10~ 3 X(Mg). The mass of synthesized 
stable nickel, M( 58 Ni)+ M( 60 Ni) = 0.006 M , is from Nomoto et al. (1991). The abundances 
in the hydrogen-rich component in SN 1987A are based on modeling of the circumstellar ring 
by Lundqvist & Fransson (1996), with C:N:0 = 1:5:4, by number. For sodium, magnesium, 
argon, and calcium we use 0.4 times the cosmic abundance in the hydrogen-rich zones, and LMC 
abundances of Fe (Russell & Dopita, 1992). The stable nickel abundance is determined from the 
Ni/Fe ratio in LMC from Russell & Dopita (1992). For stable cobalt, the solar abundance ratio, 
Co/Fe, from Anders & Grevesse (1989) is used. The explosion models only include synthesized 
elements. Therefore, if the abundance given by the explosive model is smaller than the LMC 
abundance, we estimate it from the mass fraction abundance in the hydrogen-rich component. In 
the Fe - He and Si - S regions cobalt is radioactive and its abundance therefore changes with time. 

The main physical differences in the 10H and 11E1 models are the different treatments of 
convection in the pre-supernova, and different 12 C(a,7) 16 rates. Woosley & Weaver (1986) 
employ the Ledoux criterion, in combination with semi-convection and over-shooting, while 
Nomoto & Hashimoto (1988) use the Schwarzschild criterion. 

The He - C regions in the 10H and 11E1 models are similar both in mass and in abundances. 
The oxygen-rich zones, however, differ considerably between the two models. First, the mass of 
the outer O - C region is much larger in the 10H model, compared to the 11E1 model, ~ 0.60 M Q 
and ~ 0.10 M , respectively. Secondly, in the 10H model the remainder of the oxygen zone 
has silicon and sulphur as the most abundant elements, with abundances ~ 0.065 and ~ 0.039, 
respectively. The oxygen zone in the 11E1 model is, on the other hand, dominated by neon and 
magnesium with abundances of ~ 0.11 and ~ 0.071. Only the inner part contains large fractions 
of silicon and sulphur. The mass of the O - Si - S zone is in the 10H model ~ 1.2 M , while in the 
11E1 model the O - Ne - Mg contains ~ 1.8 M©. The composition of the Si - S zone is similar in 
the two models, but the mass is in the 10H model ~ 0.30 M Q , while it is in the 11E1 model only 
~ 0.10 M.q. A case of special importance is the mixing of calcium into the oxygen zone at the time 
of oxygen burning in the 10H model. In the 11E1 model this mixing is nearly absent. The calcium 
abundance is therefore much larger in the 10H model than in the 11E1 model. This has important 
consequences for the evolution of the observed Call lines, as will be further discussed in Paper II. 

The velocity distribution and filling factors of the different components require a special 
discussion. Herant & Benz (1992) find that the central 'nickel bubble' extends to ~ 1600 km s , 
that the oxygen shell is located between 600 — 2500 km s _1 , and the helium shell between 
700 — 2600 km s _1 . The hydrogen envelope penetrates in these simulations down to ~ 600 km s~ x . 
The exact numbers are sensitive to both the progenitor model and the nature of the perturbations 
put into the simulations, and the results should therefore only be taken as indicative. However, 
a large filling factor of the iron component, coupled with similar low filling factors of both the 
oxygen, helium and hydrogen components inside of 1700 - 2000 km s _1 seem to be a generic 
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feature. In the intermediate region between 2000 - 3000 km s _1 both the hydrogen, helium and 
oxygen components have comparable filling factors. Outside of ~ 3000 km s _1 the dynamics is 
less affected, and the one-dimensional models may be fairly realistic. The quantitative properties 
of the hydrogen envelope, however, depend on the progenitor model. The line profiles can here set 
interesting constraints, as will be discussed in Paper II. 

In addition to the hydrodynamical simulations, we use, as far as possible, the observational 
information available. Li, McCray, & Sunyaev (1993) have studied the evolution of the Fe - Co 
- Ni clumps in SN 1987A. As long as the clumps are opaque to the gamma-rays they absorb 
most of the radioactive energy, causing the blobs to expand. The result is a 'nickel bubble', 
encompassing much of the volume of the core, with a low density compared to the rest of the core. 
Li et al. estimated the filling factor for these clumps to ^ 0.30, with a favored value of ~ 0.5. 
Their filling factors can not be directly compared to ours as they assumed a maximum velocity of 
2500 km s _1 for their Fe - Co - Ni clumps, while we assume a maximum value of 2000 km s . In 
our calculations we use fp e = 0.55 for the Fe - He component. This agrees with the conclusions of 
Basko (1994). 

We use a filling factor of /si-s = 9 x 10 -3 for the Si - S component in the 10H model, and 
/si_s = 3 x 10~ 3 in the 11E1 model. These filling factors give the same density in the Si - S as in 
the O - Si - S/O - Ne - Mg regions, respectively. 

Spyromilio & Pinto (1991), as well as Li & McCray (1992), find from the observed 
[OI] A A 6300, 6364 doublet ratio a local oxygen density of No(t) ~ 6.2 x 10 10 i^cm -3 . Assuming 
an oxygen mass Mo ~ 1.3 M Q , consistent with the values found from explosion models, they 
find a filling factor of the oxygen clumps of /o ~ 0.1. Assuming the same density in the 
oxygen-rich regions we find /o-Si-S = 0.06, /o-C = 0-03 for the 10H model, and /o-Ne-Mg = 0.09, 
/ _ c = 5.0 x 10~ 3 for the 11E1 model. 

For the helium-rich zone we use /n e = 0.15 within the core, and the filling factor used for the 
hydrogen within the core is /h = 0.20. These estimates are based on the simulations of Herant &: 
Benz (1992). 

For the hydrogen envelope we have tested the density structure of the Shigeyama & Nomoto 
(1990) 14E1 envelope model, as well as our own parameterized density structure, based on the line 
profiles (Paper II). The latter model is given by 

p = 9.1 X 10- 16 ( ) ( - -) g cm" 3 (38) 

r V 500 days/ \ 2000 km s" 1 / 

at V > V cove . This density is flatter than both the 11E1 and 14E1 models. 

Because the deposition, and therefore the different line luminosities, depend directly on the 
individual gamma-ray optical depths, we give these at 500 days in Table |||. For other epochs these 
scale as t~ 2 . Note that positrons add an additional amount of energy to the Fe - He region, which 
dominates the energy input to this region at late epochs. 
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3. RESULTS 

3.1. Temperature Evolution 

In Figure || the temperature evolution for the different components in the core are shown for 
the 10H and 11E1 model, respectively. 

For the metal-rich shells, Figure [2] shows an initial slow decline, followed by a steep decrease 
of the temperature during a few hundred days, after which the temperature levels off again. 
The rapid decline in temperature is due to a thermal instability, and has been referred to as the 
IR-catastrophe (Axelrod 1980, Fransson & Chevalier 1989, Fransson 1994). The instability first 
occurs in the Fe - He zone, followed by the Si - S and the O-Si-S/O-Ne- Mg zones, and 
finally by the O - C zone. 

This instability is not seen in the helium, or hydrogen-rich regions. The reason for this is 
that at the time when fine structure cooling takes over in the oxygen-rich gas, at ~ 1000 days, the 
cooling rate is a factor of ~ 4 lower in the helium-rich than in the oxygen-rich region. For hydrogen 
the corresponding factor is ~ 75. At the time when fine structure cooling in the helium and 
hydrogen-rich regions becomes important, and could cause an instability, adiabatic cooling already 
dominates. The transition to the adiabatic phase occurs at ~ 700 - 800 days in the hydrogen core, 
and at ~ 900 - 1000 days in the helium-rich core. When adiabatic cooling dominates T oc t~ 2 . 

The temperature in the time-dependent case is higher than in steady state. The reason is that 
energy from previous stages, when heating, and therefore also radiative cooling were important, is 
stored in the gas, and is then gradually lost in adiabatic expansion. In steady state a decrease in 
the heating is directly balanced by cooling. The adiabat which a given gas element follows is in 
the time-dependent case set by the epoch when radiative cooling is equal to the adiabatic. 

Because the contribution from the different zones to a given line reflects the temperature of 
this zone, we now comment on the detailed evolution of each of these zones. To guide in this 
discussion we show in Figures y and || the contributions to the fractional cooling from some of 
the most important transitions in the H-core, the He-core, the O - Si - S, and the Fe - He zones. 
Note that these figures give the collisional cooling of the gas (see equation |3|) , which is in general 
not the same as the radiative losses. In particular, transitions which are completely radiatively 
forbidden may be strong collisional coolants. An example is the a 6 D 5 / 2 — a 6 Dg/2 A 14.98 //m 
transition in the ground multiplet of Fe II, which contributes a comparable amount to the cooling 
as the observationally important a 6 L> 7 / 2 — a e Dg/ 2 A 25.99 /im transition. The excitations in the 
A 14.98 fim transition mainly decay by emitting line emission at A 25.99/iin and A 35.35/im. We 
also show in Figure |5| the relative importance of adiabatic cooling. 

In the Fe - He rich shells heating is dominated by non-thermal heating of the free electrons 
at all times. A large fraction of this is due to the deposition of positrons. At 200 days this is ~ 45 
% of the non-thermal heating, while at 1500 days the corresponding number is ~ 97 %. At most 
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20 % (but most of the time less than 10 %) of the heating is due to photoionization. The cooling 
is at all times dominated by line cooling from Fe II. The IR-catastrophe sets in at ~ 500 days, 
when the temperature is 2700 K, and is directly reflected in the relative fluxes from the iron core 
(Fig. |6|). As the temperature falls, there is a gradual progression of cooling from UV, to optical 
lines like A 7155, to near-IR, e.g., A A 1.26, 1.64 /zm, to far-IR fine structure lines. At t ^ 500 
days, Fe II A A 17.94, 25.99 /im are the most important, and finally, at t ^ 750 days, only the A 
25.99 //m line contributes to the cooling. Once the IR-catastrophe has occurred, the temperature 
is nearly constant at ~ 140 K. Because of the efficient fine-structure cooling, adiabatic cooling is 
not important before 2000 days. 

From the observations of the [Fe II] A 25.99 /im line we can get an observational estimate of 
the temperature in the iron region, assuming that the A 25.99 i/m line is the dominant coolant. 
This only applies to the iron-rich gas and not to e.g. the hydrogen- or helium-rich components. In 
the iron-rich gas the optical depth of the A 25.99^m line is 



t 26m = 28.5 /p, 1 

t 



M Fc \ / V 



500 days 



0.07 M y V2000 kms- 1 
" 2 (l- e -554K/T ) _ (39) 



With a filling factor / Fc « 0.5 (Li, McCray & Sunyaev 1993), and a temperature £ 200 K (Fig. |), 
the A 25.99 /im line is likely to be thick to at least ~ 2000 days. The critical density of the A 
25.99 fxvn line is ~ 2 x 10 4 cm -3 , so the a 6 D 9 / 2 and a 6 D 7 / 2 levels should be in LTE, and therefore 
radiate at the blackbody rate, 

dE 8irhc 



dtdV A 4 t (e hc / XkT - 1) 

2.55 x 10~ 12 (e 554 K ' T - I)' 1 

erg s _1 cm~ 3 . (40) 



t ~ l 



,500 days 

(e.g., Fransson 1994). The heating rate in the iron-rich gas is 
dE e H Oy p 

T „ 2 l K 7,56 -^7,56 + K 7,57 -^7,57j 

z ±vrji, core 

2 



dtdV 



0.75 x 10- 9 /.T 1 D~ t (— ^— ) (41) 

JFc ' V0.07 M / 
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2000 kms- 1 / \500 days, 
g-t/lU.3- (1 + L2 x 1Q -3 gt/155.6") erg 8 -l cm -3_ 

We have taken a fraction en ~ 0.5 going into heating (KF92), and assumed a 57 Co mass of 
3.3 x 10~ 3 Mq. At t ^ 900 days the A 25.99 /im line dominates the cooling (Fig. y), and we can 
estimate the temperature in the gas by balancing the two expressions above. Setting L> 7 ~ 2.5 
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(KF92) we obtain 



T e « 554 {ln[l + 1.4 x 10 -d — -^- / Fe 

1 L V0.07 M ' 
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J}- 1 K. (42) 



(1 + 1.2 x 10- 3 e */155.6^ 

At 900 days we find for f Fc = 0.5 that T e ra 191 K, at 1000 days T e « 140 K, at 1200 days T e « 97 
K, and at 1500 days T e ~ 73 K. These temperatures are about a factor of two less than those 
in Figure g. In spite of this, equation ( f42|) provides a good illustration of the sensitivity of the 
temperature to the parameters involved. It also reproduces the leveling off of the temperature in 
the Fe - He core at t ;> 1000 days (Fig. g). 

The total luminosity in the A 25.99 ^m line from the Fe - He core is 

327T 2 h C /FeKore & 
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The observed luminosity at 636 days was ~ 6.0 x 10 36 erg s _1 (Colgan et al. 1994). Multiplying 
this by a factor 1.7 to account for dust absorption, and assuming f Fe ~ 0.5 gives a temperature of 
1270 K at this epoch. However, in Paper II we will show that there is likely to be a considerable 
contribution to the flux of the A 25.99 /mi line from the hydrogen component. This temperature 
should therefore be considered as an upper limit. In our numerical model we find T e ~ 1520 K at 
636 days. Given the uncertainties in the filling factor and the dust absorption this is probably 
consistent with this temperature. 

The temperature in the iron zone has also been calculated by Li, McCray & Sunyaev 
(1993). In their calculations they consider clumps consisting of only iron, cobalt, and nickel. In 
nucleosynthesis models there is a significant fraction of helium (~ 50% by number), microscopically 
mixed with the iron, as a result of the alpha-rich freeze-out. In our calculations we find that, even 
for X(He) ~ 0.5 — 0.7, helium has a negligible effect on the conditions in the iron core. The reason 
is that only ~ 10 % of the energy is deposited in helium, and ~ 90 % in the iron peak elements. 
The temperature we find agrees well with that found by Li &; McCray. 

In the Si - S shells the temperature evolution differs somewhat between the two models. 
The instability sets in earlier in the 11E1 model, at ~ 800 days, compared to ~ 900 days for the 
10H model. The temperature at t ^ 1000 days is higher by a factor of up to ~ 1.4 for the 11E1 
model, because of a higher positron deposition in this model. Non-thermal heating dominates at 
all epochs, but photoionization of mainly Cal, Si I, and SI contribute up to as much as 40 % of 
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the heating. The positron contribution to the non-thermal heating is 10 - 20 % at 200 days, and 
more than 70 - 90 % at 2000 days. Important coolants in these shells are Si I, SI, Call, Fe I - II, 
and Ni I - II. 

Also in the O - Si - S shells in the 10H model and the O - Ne - Mg shells in the 11E1 model 
the temperature evolution differs somewhat. The instability sets in somewhat earlier in the 11E1 
model, at ~ 800 days compared to ~ 1000 days for the 10H model. We find that the temperature 
is slightly lower in the O - Si - S zones in 10H, compared to that in the O - Ne - Mg zones of 
11E1. At 2000 days the temperature is ~ 100 K in both models. 

Calculations of the temperature in the oxygen core, with and without CO, have been made 
by Liu & Dalgarno (1995). In their CO free case they find a temperature of ~ 2200 K at 800 days, 
which is considerably higher than ours. The reason is likely to be found in the differences between 
our models in the composition and in the trace elements included. Li & McCray (1992) estimate 
the oxygen temperature from observations of the [OI] A A 6300, 6364 emission doublet to T ~ 
4100 K at 200 days, and T ~ 2900 K at 500 days, in fair agreement with what we find. 

The temperature within the O - C region decreases from ~ 4000 - 5500 K at 200 days, to ~ 
2700 K at 600 days, to ~ 1500 K at 1000 days, and to ~ 110 K at 2000 days. However, the O - 
C zone is favorable to CO formation. Because we neglect the very efficient cooling by CO, our 
temperatures are probably over-estimated by a large factor. Liu & Dalgarno (1995) find that for 
the CO emitting regions the temperature is constant at ~ 1800 K for the first year, after which it 
decreases to ~ 700 K at 800 days. We will discuss the implications of this for the line emission in 
Paper II. 

In both of the oxygen-rich zones the temperature decreases nearly adiabatically once the 
IR-catastrophe has occurred, in contrast to the Fe - He region (Fig. |B|). The IR-catastrophe occurs 
slightly earlier in the O-Si-S/O- Ne - Mg regions compared to the O - C region. However, 
inclusion of CO-cooling will affect this result. 

For the helium shells within the core, line cooling is of the same order as adiabatic cooling 
after ~ 1000 days, resulting in a faster decline in the temperature than the adiabatic T oc t~ 2 
cooling. Outside of the core, adiabatic cooling becomes important at a somewhat earlier time at 
~ 800 days. 

In the hydrogen shells within the core, adiabatic cooling dominates line cooling after ~ 800 
days (Fig. |5|). Figure shows the temperatures for the different shells in the hydrogen envelope. 
The envelope density is in this, and in all other models, given by equation (pq), unless otherwise 
stated. The temperature in the envelope is at first constant, or even decreasing, up to ~ 250 days. 
This decrease is sensitive to the exact gamma-ray transfer, and therefore model dependent. After 
~ 250 days the temperature decreases nearly adiabatically, with the highest temperature closest 
to the core. Adiabaticity sets in first in the high velocity regions with the lowest density. 
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3.2. Ionization 

In Figure we show the evolution of the ionization for the same shells in the 11E1 model for 
which we previously showed the temperature and cooling. The most interesting results are the 
gradual increase in the ionization as we go to more and more advanced burning stages. This is 
mainly a result of the decrease in the number of ions for constant density as we go to the heavier 
elements. Later than ~ 1000 days the ionization in the iron core becomes nearly constant with 
time, in contrast to the oxygen-rich zones. 

Freeze-out of the ionization was discussed in connection to the bolometric light curve in 
Fransson & Kozma (1993). Because of the decreasing density the recombination time scale may 
become longer than the radioactive decay time scale. The evolution of the ionization differs 
between the different composition regions. The freeze-out effects sets in for all compositions at 
~ 800 — 900 days. However, the degree of freeze-out differs between the composition zones, and 
to a larger extent, between the core and the envelope regions. The ratio of emitted energy, L em , 
to deposited energy, L& cp , is a useful measure of the degree of freeze-out in the different regions 
(Fransson & Kozma 1993). Within the metal-rich regions in the core L em /Ld G p ~ 1-5 is reached. 
The corresponding number for the hydrogen-rich regions within the core is ~ 1.6 — 2.0. In the 
envelope, freeze-out sets in somewhat earlier (around 750 days for the regions furthest out) and a 
maximum of L em /Ld ep ~ 2.5 is reached. To see the effect of the freeze-out on the ionization we 
show in Figure [9| the electron fraction, x e , as a function of velocity for the hydrogen-rich shells at 
three epochs, 300, 800 and 2000 days. The lower ionization at V <£ 2000 km s _1 is caused by 
the higher density assumed in the core region, compared to the envelope (dotted line in Fig. ^). 
Because of the lower density, freeze-out effects are increasingly pronounced in the outer regions 
of the envelope. This is illustrated by the fact that the decrease in ionization is smaller in these 
regions than in the inner region of the envelope and within the core. 

At 1000 days, the electron fraction and temperature in the envelope are both a factor of 
2 — 4 higher in the time-dependent case, compared to steady state. In the core, the increase 
in temperature and electron fraction is ^ 1.4 for the hydrogen and helium regions, while the 
difference in the metal-rich regions is small. 



4. EFFECTS OF UV-FIELD AND CHARGE TRANSFER 

The UV-radiation field determines the ionization balance especially of ions with low ionization 
potentials and with small abundances, like Na I, Mg I, Si I and Fe I. Our model probably gives 
a good representation of the emissivity in the different lines and continua. As has been discussed 
by several authors (e.g., Xu & McCray 1991, Fransson 1994, Li & McCray 1996), this radiation 
is severely affected by scattering by the many metal resonance lines present. In addition, in 
the hydrogen-rich regions scattering by H2 may be important (Culhane & McCray 1995). The 
inclusion of these effects is beyond the scope of this paper. In Paper II we comment on their 



23 



importance. 

Besides the UV-field, charge transfer reactions are probably the most uncertain point in our 
calculations. Reaction rates involving hydrogen and helium are reasonably well known, because 
of their importance in the interstellar medium. Reactions between metals are in most cases 
uncertain, and in many cases completely unknown. In the cases where rates are available these 
are usually motivated by their importance in atmospheric physics. Unfortunately, extrapolations 
from ions with similar atomic structure are usually impossible. In some cases one can make an 
educated guess whether a certain reaction will be fast, if a near resonance between the energy 
levels of the different ions is present. This is e.g. the case for O II + Ca I — ► O I + Ca II, 
where the 5p 2 Pi/2 level in Ca II has an energy which together with the ionization potential of 
Ca I is nearly equal to the O I ionization potential. Consequently, the reaction rate is high, 
3.8 x lCP 9 (T e /10 3 K) 1 ' 2 cm 3 s _1 (Rutherford & Vroom 1972). There are also cases where this 
rule of thumb fails, and it is therefore of interest to estimate the influence of these uncertainties 
for our calculations. 

An illustration of the sensitivity to charge transfer can be obtained by a simple model with 
only two elements. Of particular interest for the discussion of the O I recombination lines in 
Paper II is charge transfer of Si I and O II in the O - Si - S component. The Si I + O II — > Si 
II + O I reaction can, like the O II + Ca I reaction, be expected to be fairly rapid, based on the 
near resonance with the 3p 4 P state in Si II, with an energy difference of only 0.12 eV from the 
difference in ionization potentials of O I and Si I. Each charge transfer reaction will then give rise 
to a photon in the Si II] A A 2334 — 2350 multiplet, contributing to the UV-field. 

In steady state the ionization balance is determined by (see eq. |10t| ) 

Toi^OI = re SiI^OII?OII+SiI + a Oll n e noil 

(r Si i + "-onCon+Sii) n Sii = asnin e nsiii. (44) 

where Toi is the non-thermal ionization rate of O I, and similarly for Si I (eq. |l(|). In addition, 
we have the number and charge conservation equations 



™OI + "-on = 


= n 


™SiI + "-Sill = 


= n Si 


nsm + n ii = 


= n e 



(45) 



These equations can easily be solved for a given epoch, and in Figure 10 we show the ionic 
abundances at 800 days for the case AT(Si) = 1.6 x 10~ 2 and X(0) = 0.98 as a function of 
the charge transfer rate between Si I + O II. These abundances are similar to the 11E1 O - 
Ne - Mg region. Although the total electron fraction only varies by ~ 25 % over the range 
£on+Sii = 10~ 15 — 10~ 8 cm 3 s _1 , the individual states of ionization vary considerably. For 
£on+Sii <J 10 -11 cm 3 s , basically all O II has disappeared, while the fraction of Si II increases. 
This behavior is easy to understand, since when £on+Sii n Sii/ a on n e <; 1 the O II fraction 
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decreases proportional to this ratio. The fraction of Si II increases correspondingly. The total 
electron fraction changes by a factor aon/^Sin- In this case aon/^Siii ~ 1) an d only a minor 
change occurs. If nsi <J n e the neutral fraction of silicon will decrease substantially. If instead 
n Si 3> n e, only a minor change in the fraction of Si I will occur. This illustrates a general 
characteristic of charge transfer. Only the minor constituents, including ions of abundant elements, 
of a given zone are affected be this process, in this case O II, Si I and Si II, but not O I. This 
agrees with the result of Swartz (1994), where he finds that if the helium abundance is larger than 
the metal abundance the charge transfer reactions do not greatly affect the ionization balance and 
level populations of helium. 

We have tested this simple model by assuming an O II + Si I charge transfer rate of 
10~ 9 cm 3 s _1 in our full code. This leads to a large decrease of the O II fraction, as expected. 
There is, however, little change in the Si I or Si II fractions. The reason for this is that another 
charge transfer reaction, Mg I + Si II — > Mg II + Si I has a large rate. Therefore, Si II only acts 
as a mediator of the charge transfer from O II to Mg II. This illustrates the gradual transfer of 
the ionization to the low ionization potential elements, as well as the general importance of charge 
transfer for the low abundance ions (including ionized states of abundant elements, e.g. O II). 



5. DUST COOLING 

The presence of dust in SN 1987A can have important consequences for the line emission. 
Besides the effect of absorption, which will be discussed in Paper II, dust cooling may also 
influence the thermal evolution of the regions where dust has formed. Although a self-consistent 
discussion of this involves a solution of the equations describing the nucleation, which is beyond 
the scope of this paper, we can see the qualitative effects of the dust by a simple model, relying on 
previous studies of the dust formation in SN 1987A by Kozasa, Hasegawa & Nomoto (1989, 1991). 

Kozasa et al. show that the most likely monomers to form are corundum, AI2 O3, and enstatit, 
Mg Si O3. AI2 O3 forms first at a temperature of ~ 1600 K, while Mg Si O3 forms at a somewhat 
lower temperature, ~ 1500 K. While Kozasa et al. mainly considered microscopically mixed cases, 
we think this is unlikely, as already argued. Therefore, there is little overlap of regions with a 
substantial abundance of iron and oxygen, and the formation of Fe3 O4 is likely to be inefficient. 
Here we consider mainly the formation of AI2 O3 and Mg Si O3 in the O - Ne - Mg zone, where 
both magnesium, silicon, and aluminum are abundant, together with oxygen. 

To model the formation of these grains we simply assume that they form at the condensation 
temperatures above. Kozasa et al. (1991) find a radius for the AI2 O3 grains of ~ 10 A, while 
the Mg Si O3 grains have a radius of ~ 70 A. The cooling rate per dust grain is, assuming that 
Tdust < T e , (Dwek & Arendt 1992) 

A dust = 1.26 x 10- 27 [ " + (^) 1 / 2 n e ] a 2 T^ erg s"\ (46) 

A I flip 



-25- 

where a is the grain radius in A. The first term represents cooling by atoms and the second that 
by the free electrons. 

The total amount of dust formed is uncertain. Here we assume, consistent with Kozasa et al. , 
that all aluminum in the O - Ne - Mg zone goes into AI2 O3 and all silicon into Mg Si O3 grains. In 
the 11E1 model, X(Mg) « 7.1 x 10~ 2 , and A(Si) « 1.6 x 10~ 2 . We take the aluminum abundance 
from the 11E1 model to be 8.3 x 10~ 3 by mass. The total mass of the O - Ne - Mg zone is 1.8 M , 
giving M(A1) = 1.5 x 10~ 2 M Q , and the mass of Al 2 O3 is then 2.8 x 1CT 2 M and the mass of 
Mg Si 3 ~ 0.17 M . With a size of 1.7 A and 2.06 A per key species (Kozasa et al. 1989), the 
total number of grains in the two cases are therefore 1.8 x 10 51 (M(A1)/1.5 x 10~ 2 M Q ) (a/10 A)~ 3 
and 5.2 x 10 49 (M(Si)/5 x 10~ 2 M ) (a/70 A)" 3 , respectively. The total, optically thin cooling 
rate for AI2 O3 is then 

A o _ llxl0 -i3 (M(0)\ ( M(A1) \ 
A dust - 1.1x10 \jy^) { 1 , 5xW -2 Mq ) 

V \ - 6 / i \ - 6 / f \ - 2 

''core \ / b \ I JO 



and for Mg Si 3 
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„ \ -1 
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dust 
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As the temperature falls the cooling will eventually exceed the black-body limit, at which time the 
temperature stabilizes. 

Assuming that dust formation occurs when the temperature has decreased to the condensation 
temperatures for the two species, the AI2 O3 formation takes place at ~ 800 days. At the time 
of dust formation in our model, immediately before the temperature drop, the ratio of the dust 
cooling and line cooling is ~ 50. Therefore, as soon as the dust forms it immediately cools the 
gas to a very low temperature, at which the optically thick dust cooling balances the heat input. 
Because of the dramatic cooling, the Mg Si O3 formation occurs at the same epoch. 

A dust formation epoch at ~ 800 days is later than the observations indicate. Possible reasons 
for this discrepancy could either be extra cooling in the oxygen-rich gas due to e.g. molecular 
cooling by SiO, or the simplified assumption of instantaneous dust formation at the temperatures 
above. Other possible, and perhaps the most likely, sites for dust formation are in the O - C 
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zone and in the iron-rich core. In the O - C region CO cooling could decrease the temperature 
enough for graphite to form. This scenario requires that not all carbon is locked up into CO. 
This is supported by the calculations done by Liu & Dalgarno (1995). They find a CO mass of 
(1 - 2) x 10~ 3 M Q , while the total carbon mass in the O - C zone is (1 - 10) x 1CT 2 M Q . A 
second possible site is the iron core, where the temperature is <^ 1500 K at ~ 600 days. Clumping 
may further decrease the temperature and thus the dust condensation epoch. The dust formation 
is likely to occur at different epochs at different positions in the ejecta, explaining the observed 
gradual formation of the dust. 

In the optically thick limit we can calculate the dust temperature by balancing the gamma-ray 
energy absorbed in the O - Ne - Mg zone by the dust cooling to obtain 



V \~ 1/2 / t 

1770 ' 



2000kms-V V 500 days 
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We have here scaled the gamma-ray optical depth of the O - Ne - Mg zone by the value in our 
models at 500 days (Table ||). At 600 days we obtain Td us t = 380 K, and at 800 days 183 K. These 
values are in reasonable agreement with those determined by Wooden et al. (1993) and Colgan et 
al. (1994). 

Although certainly too simplified, this discussion is indicative of the importance of the 
additional cooling where dust is formed. The cooling rate we estimate above is directly proportional 
to the amount of dust assumed to be formed, as long as the gas is optically thin. Because of the 
nearly full conversion of the key elements into dust, this should probably be considered as an 
upper limit. A serious discussion of the dust cooling requires a more physical model of the dust 
nucleation, along the lines of Kozasa et al. 

A consequence of this scenario is a depletion of the elements in the oxygen region involved 
in the dust condensation. This has little effect on oxygen. Both aluminum and silicon are by 
assumption completely depleted from the gas phase in this region, and magnesium is decreased by 
~ 20% in the 11E1 model. Observational evidence for Si depletion, based on the [Si I] lines, has 
been prestented by Lucy et al. (1991). It is, however, not trivial to separate the effects of dust 
depletion and decrease in temperature on the strength of the Si emission. 



6. CONCLUSIONS 

The temperature and ionization of the supernova ejecta determine the evolution of the 
observed emission lines. In this paper we have made a considerable effort to calculate these as 
realistically as possible. Previous calculations have made various simplifications, which limit their 
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applicability. Fransson & Chevalier (1987, 1989) neglected time-dependent effects, as well as 
the hydrogen-rich component. Li, McCray & Sunyaev (1993) and Liu & Dalgarno (1995) only 
considered the iron and oxygen-rich zones, respectively, and also neglected time-dependent effects. 
The only work comparable to ours is that by de Kool, Li &; McCray (1997), which we became 
aware of only after this work was completed. In as far as we have been able to compare, we agree 
very well with their results in terms of both temperature and ionization. This comparison of 
completely independent calculations adds to the confidence of our results, within the assumptions 
of the model. 

Except as input to the calculations of the line emission, our most important results concerns 
the IR-catastrophe at 600 - 1000 days in the metal zones. Other results of special interest concern 
adiabatic cooling later than 500 - 800 days in especially the hydrogen and helium regions, and the 
freeze-out of the ionization in the hydrogen envelope. Although some of these effects have been 
noted previously (Fransson & Chevalier 1987, 1989; Fransson & Kozma 1993), they are in this 
paper synthesized into a coherent picture. In Paper II we show how the temperature evolution is 
reflected in the line emission. In addition, both molecule and dust formation are sensitive to the 
temperature. In our models dust is formed too late in the O - Ne - Mg zones. 

In spite of our effort to include physics as realistically as possible, there are still some major 
deficiencies, in particular our treatment of the UV-field, the neglect of molecular and dust cooling, 
and uncertainties in the charge transfer rates. Molecular cooling can have important consequences, 
as demonstrated by Liu & Dalgarno (1995), but is probably limited to the O - C zone, and the 
interface between the oxygen and silicon regions. As we will discuss in Paper II, the thermal line 
emission from the other regions would otherwise also be completely quenched. As our estimate 
of the dust cooling shows, a similar effect would occur where the dust forms. In this case, 
however, dust formation takes place where there is a strong drop in the temperature, because of 
the IR-catastrophe. The emission from the metal-rich zones is then dominated by non-thermal 
excitation, and the line emission is therefore insensitive to the temperature. These areas clearly 
deserve much further study. 
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A. APPENDIX 

A.l. Charge Transfer 

Charge transfer reactions are important in determining the ionization balance for the trace 
elements. Unfortunately, rates are rather uncertain. The reactions we include in our calculations 
are given in Table || together with references for the rates. Most charge transfer rates with 
hydrogen are from Kingdon & Ferland (1996), who give analytical fits to all the rates. This 
includes a compilation of both already available rates from different sources (see references 
therein), and newly calculated rates. 

The data in Prasad & Huntress (1980) and Arnaud & Rothenflug (1985) are also compilations 
from various sources. The rates from Kimura, Lane, Dalgarno, & Dixson (1993) (H II + He I — * 
H I + He II), and Kimura et al. 1993 (He II + C I -> He I + C II) are low. For the reaction H II 
+ Ca I — ► H I + Ca II we did not find any published data. As an estimate we use the same rate 
as for O II + Ca I -> O I + Ca II. 

Swartz (1994) estimates charge transfer reaction rates between helium and metals, using the 
Landau-Zener and modified Demkov approximations. He also estimates the rates for the charge 
transfer between excited states in He I (2s 1 S, 2p 1 P, 2s 3 S, and 2p 3 P) and metals. The accuracy 
of these are probably not better than a factor of two to five. We therefore include the rates from 
Swartz separately, to examine their importance. 

The rate for the Penning process, 7.5 x 10" 10 (T/300K) 1 / 2 cm 3 s- 1 , is from Bell (1970). 



A. 2. Model Atoms 

A. 2.1. HI 

For H I we use a model atom consisting of 30 levels. For the 5 lowest n-states all Z-states are 
included. The higher n-states, up to n = 20, are treated as single levels. 

We use expressions given by Brocklehurst (1971) for the transition probabilities between 
different nZ-states (A n i n >i/). Only transitions with Al = ±1 are assumed to have A-values greater 
than zero. The one exception is the two-photon transition, 2s — Is, with a transition probability 
of 8.23 s . When calculating transition probabilities between an upper n-state and a lower 
nZ-states, as well as between two n-states, the following averages are used, 

n-l 
An,n'V = z2 ~Anl,n'l'i (Al) 

and 



i=o dn 



A n , n ' =22 z2 Anl,n'V = /2 An,n'V (A2) 
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where g\ = 2(2/ + 1) is the statistical weight of the Estate, and g n = 2n 2 is the statistical weight 
of the n-state. In thermodynamic equilibrium, when the /-states are populated according to their 
statistical weights (x n i = x n ®-), we have 

x n-^-n,n'V = / , x nl-^-nl,n'V (■"■«*) 

I 

and 

x n-^n,n' = / t x nl-^-nl,n'l' i (.-"-4) 

1,1' 

motivating the expressions above. 

Recombination coefficients to the n/-states are calculated using expressions in Brocklehurst 
(1971). The total recombination coefficient for hydrogen is given by Seaton (1959). These rates 
agree with the total recombination coefficients by Verner & Ferland (1996). In order to get the 
correct total recombination coefficient we increase the recombination coefficient to the highest 
(n = 20) n-state, so that the sum of the recombination coefficients to all the levels equals the total 
recombination coefficient. This treatment is not correct for the highest levels, but should be a 
good approximation for the lower levels. 

Collisions with both electrons and protons between /-states for a given n-level are included for 
the 5 lowest n-levels (Brocklehurst 1971, Pengelly & Seaton 1964). The total collisional excitation 
rates between different n- levels (C nn i) are from Johnson (1972). To divide the total rates into 
rates between different /-states we use the Bethe-approximation (Hummer & Storey 1987). For all 
An > and A/ = ±1 we use, 

(^nl,n'V = "H C nn ' (A5J 

J\ r , 



'mi 



and 



This approximation gives, 



Cn,n'l' — ~~T (-'nn'- (A6) 



^nrv 



x n^n,n' — 2_^ X nl^- / nl,n'V i {"■') 



l.V 

if the different /-states are assumed to be populated according to statistical weights. 

Photoionization cross sections for individual nZ-states, for n > 2, are calculated using 
the expressions given by Brocklehurst (1971), and for the higher n-states by Seaton (1959). 
Photoionization cross sections for n = 1 and n = 2 are from Bethe & Salpeter (1957). 

To solve the level populations for the hydrogen atom using the Sobolev approximation 
complications appear due to overlapping lines. Transitions between two n-states but different 
/-states have the same energy, and the lines may interact. However, because of the degeneracy 
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this scattering is local, and the escape formalism may still be used. This case has been treated by 
Drake & Ulrich (1980), and we use their expressions. 

The treatment of an n-state as a single level is correct for high enough densities, in which case 
the /-states become populated according to statistical weights. When calculating the line emission 
we recalculate the level populations for hydrogen using a 210 level atom containing all /-states up 
to n = 20. This is done separately from the calculation of temperature and ionization, where we 
assume /-mixing for n > 6. 

For the hydrogen two-photon continuum we use the two-photon distribution given by 
Nussbaumer & Schmutz (1984). 



A. 2.2. He I 

We use a model atom for He I, consisting of 16 levels. The first 12 levels are Is 1 5', 2s 3 S, 
2s l S, 2p 3 P, 2p l P, 3s 3 S, 3p 3 P, 3d 3 D, As 3 S, Ap 3 P, Ad 3 D, Af 3 F. The remaining 4 levels 
are fictitious levels, into which all the triplet levels up to n = 100 have been combined. Collision 
strengths, and transition probabilities are from Almog & Netzer (1989). Verner & Ferland (1996) 
have calculated recombination coefficients, both total and to the different excited states, for the 
entire temperature range, 3 - 10 9 K, which we use. 

In our model we use the direct recombination coefficients calculated by Verner & Ferland 
(1996) to all levels except to 2s 1 S and 2p 1 P and the Ap 3 P and Af 3 F states. The reason for 
this is that the 2s 1 S and 2p l P states are the only excited singlet levels included in the model 
atom by Almog & Netzer (1989). Including only direct recombinations to these would therefore 
under-estimate the total recombination rate, including higher states. We use the fit given by 
Almog & Netzer, which is, however, uncertain at low temperatures. This fact introduces some 
uncertainty in the He I two-photon continuum at low temperatures. 

For the Ap 3 P and 4/ 3 F levels we use rates from Almog & Netzer data in the temperature 
interval 10000 to 20000 K, extrapolating this to lower temperatures by a a T -0 ' 5 . This 
is in agreement with the temperature dependence Bates (1990) finds for low temperatures. 
Recombination coefficients to the four highest, fictious levels are from Almog & Netzer in the 
temperature interval 10000 - 20000 K. The temperature dependence is corrected for lower 
temperatures, using a oc T~ 0,5 for T < 1000 K. The rates of these four highest recombination 
coefficients were normalized to achieve the correct total recombination coefficient. 

For the ground state photoionization cross sections are from Verner et al. (1996), while the 
excited state cross sections are from Koester et al. (1985). For the helium two-photon continuum, 
between the 2s *S and the Is X S states, we use the distribution given by Drake, Victor, Sz Dalgarno 
(1969). 

In the Spencer-Fano calculations we include non-thermal excitations to the following excited 
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levels: 2s 3 S, 2s 1 S, 2p 3 P, 2p 1 P, 3s 3 S, 3p 3 P, 3d 3 D. References for the cross sections are given 
in KF92. 



A. 2.3. I 

Our 01 atom consists of 13 levels, 2p A 3 P 2 , 2p 4 3 P ± , 2p A 3 P , 2p A 1 D, 2p A 1 S, 3s 5 S" , 3s 3 S , 
3p 5 P , 3p 3 P , 4s 5 S , 4s 3 S , 3d 5 D , 3d 3 D . Transition probabilites and collision strengths 
are from Bhatia & Kastner (1995). We include non-thermal excitations from the ground state to 
all excited levels. References for the cross sections used are given in KF92. For the fine structure 
transitions within the 3 P state A-values and collision strengths for O I are from Berrington (1987). 

The total radiative recombination coefficient is from Chung, Lin, & Lee (1991). Recombination 
coefficients to the individual levels are from Julienne, Davies &; Oran (1974), renormalized by a 
factor 1.38 to give a total rate equal to that of Chung et al. 

Photoionization cross sections are from Verner et al. (1996) for the ground state. For the 
excited states the data are from Dawies & Lewis (1973). 



A.2.4. Ca II 

The six-level model atom for Call consists of the 4s 2 S, 3d 2 D, 4p 2 P, and 5s 2 S levels. 
A-values are from Ali &; Kim (1988), Zeippen (1990), and Wiese et al. (1969). Collision strengths 
are from Burgess et al. (1995), except between the sublevels of 3d, which are from Shine (1973). 
Ground state photoioinization cross sections are from Verner et al. (1996). For the excited states 
data are from Shine (1973), while recombination rates are from Shull & Van Steenberg (1982). 



A. 2. 5. Fe I-IV 

Fe I-IV are treated as multilevel atoms with 121 levels for Fel, 191 levels for Fell, 110 levels 
for Felll, and 43 levels for FelV. A-values for Fel are from Kurucz & Peytremann (1975) and 
from Axelrod (1980), and collision strengths from Axelrod (1980). The accuracy of the latter 
are probably low. For Fell, A-values and collision strengths are from the Iron Project (Nahar 
1995, Zhang & Pradhan 1995), and from Garstang (1962). Additional A-values are from Kurucz 
(1981). For Felll and FelV data come from Garstang (1957) and Kurucz & Peytremann (1975), 
Berrington et al. (1991). Total recombination rates are from Shull &; Van Steenberg (1982), and 
the fractions of the radiative recombination going to the ground state are from Woods, Shull, & 
Sarazin (1981). 

For Fel and Fell the total photoionization cross sections are from Verner et al. (1996). Only 
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photoionizations from the first five levels (i.e., the ground multiplet) are included for Fel and 
Fell. The same cross section is assumed for all levels in the ground multiplet. 

It is highly time consuming to calculate the level populations for the large Fe I-IV ions. In 
order to speed up the calculations we reduce the number of levels in these, if possible. The first 
time steps of the calculation always use 'full' iron-atoms. At later epochs, the populations of 
the 'full' atoms are calculated every few days. The highest levels, which in these calculations 
together contribute less than 0.1 % of the total emission or cooling/heating, are then excluded in 
the following time steps. The error in the total emission introduced by this procedure is therefore 
never larger than ~ 0.1%. 



A. 2.6. Nil- II and Co II 

We include Ni I, Ni II, and Co II with a small number of levels in order to model the [Ni I] 
A 3.119 /on, [Ni I] A 7.505 /mi, [Ni II] A 6.634 /mi, [Ni II] A 10.68 /mi, and the [Co II] A 
10.52 /mi lines. We assume that the ionization balances for cobolt and nickel are coupled to the 
ionization balance of iron by charge transfer (Li, McCray, & Sunyaev 1993). Due to lack of good 
atomic data we do not include non-thermal excitations or recombinations in calculating the line 
emission. Therefore, lines arising from these atoms are due only to thermal, collisional excitations 
by electrons. 

Our Ni I atom has seven levels, consisting of the three lowest multiplets, a s F, a 3 D, and 
a}D. Transition probabilities between levels within the a?F and the a?D multiplets are from 
Nussbaumer & Storey (1988), and transition probabilities between a}D and a?D from Garstang 
(1964). 

The model atom for Ni II has eight levels, consisting of the three lowest terms, 2 D, 4 i ? , and 
2 F. Transition probabilities are from Nussbaumer & Storey (1982), and collision strengths are 
from Bautista & Pradhan (1996). For Co II we only include the lowest multiplet, a 3 F, consisting 
of three levels, with transition probabilities from Nussbaumer & Storey (1988). 

Unfortunately, no collision strengths are available for Ni I and Co II. Axelrod (1980) estimated 
collision strenths from the empirical relation 

ttij = Cg igj (A8) 

where gi, and g« are the statistical weights for the levels and C is a constant. The constant was 
determined by a comparison with calculations by Garstang, Robb & Rountree (1978) for Fe III, 
and found to be C = 4 x 10~ 3 for transitions with A < 10 /mi. For infrared transitions with 
A > 10 /mi he adopted a higher value of the constant, C = 2 x 10~ 2 . For the lowest 16 levels 
in Fe II we have checked this relation with the new Opacity Project collision strengths (Zhang 



& Pradhan 1995), and found C = 0.046. In order to test equation (A8), we calculated collision 
strengths for Ni II using both C = 4 x 10 -3 , and C = 0.046, and compared to the collision 
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strengths from Nussbaumer & Storey (1982). The Ni II data are best reproduced by the constant 
from Axelrod, while our larger C gives a factor of 10 too large collision strengths. Therefore the 
constant from Axelrod gives acceptable values for Ni II, while it under-estimates the collision 
strengths for Fe II by a factor of 10. 

Li, McCray, & Sunyaev (1993) proposed an empirical formula for the collision strengths, based 
on known transitions in Fe II. In their formula the collision strength depends on the transition 
probability, and the energy difference between the levels, as well as on the statistical weight of 
the upper level. We do not find any clear correlation between transition probabilities, energy 
differences and collision strengths. 

Given the result of the tests above, we feel that it is not possible to estimate the unknown 
collision strengths empirically to any degree of certainty. We therefore simply put the unknown 
collision strengths equal to 0.1 and await new atomic data. 



A. 2.7. Other Ions 

References for the collisional ionization cross sections used for calculating the non-thermal 
ionization rates, are given in KF92. For the photoionization cross sections we use the analytical 
fits given by Verner et al. (1996). 

Verner & Ferland (1996) give analytical fits to the radiative recombination coefficients in the 
temperature interval 3 K to 10 9 K for recombinations towards H-, He-, Li- and Na-like ions. We 
use their data for recombination of He III, Mg III, and Na II. Radiative recombination coefficients 
from Gould (1978), and low-temperature dielectronic recombination coefficients from Nussbaumer 
& Storey (1983, 1986) are used for recombination of C II, C III, N II, O III, Mg II, Si II, Si 
III. We use rates from Shull & Van Steenberg (1982) for the Ne II, Mg II, S III, Ar II, Ca II 
recombinations. For S II we take radiative recombination rates from Gould (1978). 

To calculate the cooling and line emission, the level populations for a number of model atoms 
are calculated, in steady state. C I, N II, O III, and S III are treated as six-level atoms, and C II, 
N I, O II, Si II, S I and S II as five-level atoms. References for the atomic data used in these model 
atoms are given in Lundqvist & Fransson (1996). In addition, for Oil we include non-thermal 
ionization of 01 to the excited 2p 3 2 D° and 2p 3 2 P° levels in Oil (KF92). 

The following transitions are calculated in the two-level approximation : CII A 1334, 
[Nell] A 12.814 /im, Nal A 2843, Mgl] A 4571. A-values and collision strengths are from Mendoza 
(1983), Fabrikant (1974), Hayes & Nussbaumer (1984), Nussbaumer & Storey (1986), and Morton 
(1991). 

For Si I we include the 3p 2 3 P, 3p 2 1 D, and 3p 2 1 S levels, plus the fine structure transitions 
of the 3 P ground state. A-values are from Mendoza (1983), and collision strengths for 3p 2 3 P to 
3p 2 1 D from Pindzola, Bhatia, & Temkin (1977). For the other transitions we estimate collision 
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strengths from C I and O I. 

For Nal A 2843, Mgl A 2852, and Mgll A A 2795,2802 we include non-thermal excitation 
with cross sections given in KF92. The cross section for Mgl] A 4571 is from Fabrikant (1974). 

Lines arising as a result of radiative recombination can be strong for abundant elements, 
unless charge transfer is more efficient. An important example is recombination to Mg I, where 
the strongest line is expected to be Mg I] A 4571. No detailed calculations of the radiative 
recombination cascade of this element exist. Dielectronic recombination is unimportant, although 
we include it. Because most recombinations to the triplet levels should pass through the 3p 3 Po 
level, we estimate the effective rate as proportional to the statistical weight of the triplets, i.e. 3/4 
of the total rate to the excited states, or 2.0 x 10~ 13 (T/10 4 K) 0,86 cm 3 s -1 . The contribution of 
dielectronic recombination to the Mgl] 4571 line is from Nussbaumer & Storey (1986). 
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Table 1. Abundances by number in model 10H 



Fe-He Si-S O-Si-S O-C He-C H 

Mass 0.07 0.30 1.20 0.60 2.00 9.00 



H 


0.00 


0.00 


0.00 


0.00 


0.00 


8.00 


x IO" 1 


He 


5.56 


x IO" 1 


1.48 


x 10~ 15 


5.64 


x 


io- 12 


6.29 x 


io- 10 


9.75 


x 


io- 1 


2.00 


x IO" 1 


C 


2.75 


x 10- 4 


3.81 


x 10~ 6 


1.68 


x 


io- 2 


1.59 x 


io- 1 


1.82 


x 


io- 2 


2.56 


x IO" 5 


N 


0.00 


0.00 


0.00 


0.00 


6.88 


x 


io- 6 


1.28 


xlO" 4 





3.95 


x 10~ 7 


9.38 


x 10~ 5 


8.53 


x 


io- 1 


8.24 x 


io- 1 


2.03 


X 


10- 3 


1.04 


x IO" 4 


Ne 


4.77 


x 10~ 7 


1.17 


x 10~ 9 


9.33 


x 


10~ 3 


1.62 x 


io- 2 


3.97 


X 


10- 3 


4.40 


x IO" 5 


Na 


7.66 


x 10~ 9 


1.08 


x 10~ 7 


9.06 


x 


10~ 5 


6.22 x 


io- 6 


1.03 


X 


io- 6 


6.40 


x IO- 7 


Mg 


9.38 


x 10~ 7 


1.32 


x 10~ 5 


1.11 


x 


io- 2 


7.62 x 


io- 4 


1.26 


X 


io- 4 


9.60 


x IO- 6 


Si 


1.10 


x 10~ 5 


6.30 


x IO" 1 


6.46 


x 


io- 2 


3.79 x 


io- 4 


1.04 


X 


io- 4 


1.36 


x IO" 5 


S 


3.42 


x 10~ 5 


2.91 


x IO" 1 


3.88 


x 


io- 2 


1.90 x 


io- 4 


5.21 


X 


10- 5 


4.48 


x IO- 6 


Ar 


5.46 


x 10~ 5 


4.06 


x 10~ 2 


5.25 


x 


10- 3 


4.60 x 


10- 5 


1.26 


X 


10- 5 


2.32 


x IO- 6 


Ca 


1.91 


xlO- 4 


2.58 


x IO" 2 


5.48 


x 


io- 4 


2.82 x 


10- 5 


1.90 


X 


io- 6 


7.20 


x IO- 7 


Fc 


4.43 


x IO" 1 


1.20 


x 10~ 2 


2.20 


x 


io- 4 


1.90 x 


io- 4 


5.30 


X 


io- 5 


2.00 


x IO" 5 


Co 




- 




- 


5.50 


x 


io- 7 


4.75 x 


io- 7 


1.30 


X 


io- 7 


6.25 


x IO" 8 


Ni 


3.56 


x IO" 2 


9.62 


x 10- 4 


1.40 


x 


10- 5 


1.20 x 


io- 5 


3.40 


X 


io- 6 


1.60 


x IO- 6 



Note. - - The first row gives the mass (in Mq) for each composition region. 



41 



Table 2. Abundances by number in model 11E1 



Fe-He 



Si-S 



O-Ne-Mg 



O-C 



He-C 



H 



Mass 



0.07 



0.10 



1.80 



0.10 



2.00 



9.00 



H 


0.00 


0.00 


0.00 


0.00 


0.00 


8.00 


x IO" 1 


He 


7.32 x IO" 1 


1.73 x IO" 11 


8.47 


X 


io- 7 


5.90 


X 


10~ 4 


9.81 


X 


io- 1 


2.00 


x IO" 1 


C 


6.02 x 10~ 4 


3.64 x 10~ 6 


8.35 


X 


10~ 3 


2.71 


X 


io- 1 


1.41 


X 


io- 2 


2.56 


x IO" 5 


N 


0.00 


0.00 


0.00 


0.00 


4.76 


X 


io- 6 


1.28 


xlO- 4 





8.70 x 10~ 7 


1.92 x 10~ 5 


7.83 


X 


io- 1 


7.15 


X 


io- 1 


1.93 


X 


10- 3 


1.04 


x IO" 4 


Ne 


3.81 x 10~ 7 


3.31 x IO" 10 


1.14 


X 


io- 1 


1.06 


X 


io- 2 


2.63 


X 


10- 3 


4.40 


x IO" 5 


Na 


1.28 x 10~ 10 


0.00 


5.18 


X 


io- 4 


5.90 


X 


io- 6 


1.70 


X 


io- 6 


6.40 


x IO- 7 


Mg 


1.02 x 10~ 6 


1.44 x 10~ 5 


7.12 


X 


io- 2 


2.72 


X 


10- 3 


2.40 


X 


10- 5 


9.60 


x IO- 6 


Si 


1.47 x 10~ 5 


5.87 x 10- 1 


1.64 


X 


io- 2 


1.30 


X 


io- 4 


3.60 


X 


10- 5 


1.36 


x IO" 5 


S 


7.36 x 10~ 5 


2.76 x IO" 1 


2.10 


X 


io- 4 


4.20 


X 


10- 5 


1.20 


X 


10- 5 


4.48 


x IO- 6 


Ar 


1.86 x 10~ 4 


4.04 x 10~ 2 


2.40 


X 


10- 5 


2.10 


X 


10- 5 


6.00 


X 


io- 6 


2.32 


x IO- 6 


Ca 


6.56 x 10~ 4 


2.64 x IO -2 


7.70 


X 


io- 6 


6.70 


X 


io- 6 


1.90 


X 


io- 6 


7.20 


x IO- 7 


Fc 


2.39 x IO" 1 


6.94 x IO" 2 


2.10 


X 


io- 4 


1.80 


X 


io- 4 


5.20 


X 


10- 5 


2.00 


x IO" 5 


Co 


- 


- 


5.25 


X 


io- 7 


4.50 


X 


io- 7 


1.30 


X 


io- 7 


6.25 


x IO" 8 


Ni 


2.02 x 10~ 2 


5.69 x 10~ 3 


1.40 


X 


10~ 5 


1.20 


X 


io- 5 


3.40 


X 


io- 6 


1.60 


x IO- 6 



Note. - - The first row gives the mass (in Mq) for each composition region. 
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Table 3. Optical depths to the gamma-rays at 
500 days. 



Zone At 7 



Fe - He 0.015 

Si - S 0.061 

O - Si - S 0.225 

O - C 0.103 

He (core) 0.090 

He (envelope) 0.030 

H (core) 0.676 

H (envelope) 0.302 

Total 1.502 
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Table 4. Charge Transfer Reactions 



Reaction 


Reference 


Reaction 


Reference 


H I + He II - 


-> H II + He I 


1 


H I + He III 


-^ H II + He II 


1 


H I + C II -► 


H II + C I 


1 


H I + C III - 


-> H II + C II 


1 


H I + N II -► 


H II + N I 


1 


H I + II -i 


. H II + I 


1 


H I + III - 


-> H II + II 


1 


H I + Mg III 


-^ H II + Mg II 


1 


H I + Si III - 


-> H II + Si II 


1 


H I + S II -► 


H II + S I 


1 


H I + S III - 


► H H + S II 


1 


H I + Fe III 


^ H II + Fe II 


1 


H I + Fe IV - 


-+ H II + Fe III 


1 


H I + Fe V - 


♦ H II + Fe IV 


1 


H II + He I - 


-> H I + He II 


4 


H II + C I -» 


■ H I + C II 


1 


H 11 + N I -► 


H I + N II 


1 


Hn + oi-i 


. H I + II 


1 


H II + Na I - 


-> H I + Na II 


5 


H II + Mg I - 


-> H I + Mg II 


1 


H II + Mg II 


-> H I + Mg III 


1 


H II + Si I - 


> H I + Si II 


2 


H II + Si II - 


-> H I + Si III 


1 


H II + S I -> 


H I + S II 


1 


H II + Ca I - 


-> H I + Ca II 


9 


H II + Fe I - 


+ H I + Fe II 


3 


H II + Fe II - 


-► H I + Fe III 


1 


He I + C III 


-^ He II + C II 


2 


He I + III 


-> He II + II 


2 


He II + C I - 


-> He I + C II 


6 


He II + C II - 


-> C III + He I 


2 


He II + Si I - 


^ He I + Si II 


3 


He II + Si II 


-► Si III + He I 


2 


He II + S II - 


-► S III + He I 


2 


C II + Na I - 


-> Na II + C I 


3 


C II + Mg I - 


-> Mg II + C I 


3 


C II + Si I -> 


Si II + C I 


3 


C II + S I -»• 


S II + C I 


3 


C II + Fe I - 


> Fe II + C I 


3 


N II + Mg I - 


-> N I + Mg II 


3 


N II + Ca I - 


-> N I + Ca II 


7 


N II + Fe I - 


+ N I + Fe II 


3 


II + Ca I - 


-> I + Ca II 


8 


II + Fe I - 


-> I + Fe II 


3 


Na I + Mg II 


^ Na II + Mg I 


3 


Na I + Si II - 


^ Na II + Si I 


3 


Na I + S II - 


♦ Na II + S I 


3 


Na I + Fe II 


^ Na II + Fe I 


3 


Mg I + Si II - 


-> Mg II + Si I 


3 


Mg I + S II - 


-> Mg II + S I 


3 


Si I + S II -» 


Si II + S I 


3 


Si II + Fe I - 


-> Si I + Fe II 


3 


S II + Fe I -> 


S I + Fe II 


3 









References. — (1) Kingdon & Ferland 1996; (2) Arnaud Sz Rothenflug 1985; (3) Prasad & Huntress 
1980; (4) Kimura, Lane, Dalgarno, & Dixson 1993; (5) Croft & Dickinson 1996; (6) Kimura et al. 
1993; (7) Swartz 1994; (8) Rutherford k Vroom (1972); (9) Estimate; the same rate as for Ca I + O 
II 
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H/He envelope 




Fig. 1. — Schematic figure showing the model structure of the ejecta, where the thickness of the 
zones correspond to the filling factors of the shells. We have for clarity omitted the O - C and Si 
- S zones, as well as the various zones in the envelope. In the white zone between the second and 
third iron zones we have an additional set of O, He, H, Fe, O, He, and H zones. The inner hydrogen 
zone and the outermost iron zone, which break the regularity, have been added to produce the 
correct line profiles. 



45 





_ I I I I I 


I I I I I I 


1 _ 






Fe-He 


- 


«, 




10H 


- 






1 1F1 




u 




£ 1000 









(0 






— 


u 






— 


CD 






— 


ft 






— 


6 






- 


V 








E-h 






- 


100 


—I I I I I I 


ill 1 1 1 1 1 1 


1 




500 1000 1500 2000 

Time (days) 



500 1000 1500 2000 

Time (days) 



"i i i i i i i i i i i i i i i i r 



O-Si-S/ 
O-Ne-Mg 




100 I — i i I i i i 




500 1000 1500 2000 

Time (days) 



100 La i I i i i i I i i i i I i i l_i 

500 1000 1500 2000 

Time (days) 



u 

3 1000 

-!-> 

«J 

u 

CD 

ft 

a 

CD 

E- 



I I I I I I I I I I I I I I I I I 



He 



100 hilin 



i i i i i i i i 



CD 

3 1000 

■<-> 

cC 

u 

CD 
ft 



CD 
E-> 




500 1000 1500 2000 

Time (days) 



100 H i I i i i i I i i i i I i i i t 

500 1000 1500 2000 

Time (days) 



Fig. 2. — The temperature evolution for a representative shell for the different composition zones 
in the core for the 10H (solid line) and 11E1 (dashed line) models. 
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Fig. 3. — Fractional cooling by the most important transitions in the H-core zone and the He-core 
zone, as a function of time. 
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Fig. 4. — Fractional cooling in the O - Ne - Mg and the Fe - He zones, as a function of time. 
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Fig. 5. — Ratio of adiabatic to radiative cooling. The hydrogen curves are from the top the 
hydrogen in the core (solid), and the hydrogen envelope at 2000 km s _1 (dotted), 2800 km s _1 
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Fig. 6. — The strongest Fe II lines from the iron core, together with the temperature in the same 
region. Note the gradual transition from optical to near-IR, to far-IR lines, as the IR-catastrophe 
sets in at ~ 600 days. 
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Fig. 7. — The temperature evolution for the different shells in the hydrogen envelope. The shell 
furthest in has the highest temperature. The velocities of the shells are ~ 2200, 2600, 3000, 3700, 
4300, 4700, 5250, 5750 km s -1 . 
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Fig. 8. — The evolution of the electron fraction for representative shells for the different composition 
zones in the core for the 11E1 model. 
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Fig. 9. — Electron fraction as a function of velocity for the hydrogen-rich regions at three different 
epochs. At 300 days steady state is a good assumption, at 800 days the freeze-out effects start 
to become important. Freeze-out effects are more pronounced in the outer envelope region (V 
> 2000 km s^ 1 ) where the density is lower. Also shown, as the dotted curve, is the number density 
for the hydrogen-rich regions at 800 days. 
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Fig. 10. — Electron and ionic fractions of O II and Si II, as a function of the Si I + O II 
+ 01 charge transfer rate. 
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